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Abstract. Ab initio formulations of the interlayer exchange coupling (IEC) be- 
tween two, in general non-collinearly aligned magnetic slabs embedded in a non- 
magnetic spacer are reviewed whereby both the spacer and the magnetic slabs as 
well as their interfaces may be either ideal or random. These formulations are based 
on the spin-polarized surface Green function technique within the tight-binding lin- 
ear muffin-tin orbital method, the Lloyd formulation of the IEC, and the coherent 
potential approximation using the vertex-cancellation theorem. We also present an 
effective method for the study of the temperature dependence of the IEC. The peri- 
ods, amplitudes, and phases are studied in terms of discrete Fourier transformations, 
the asymptotic behavior of the IEC is briefly discussed within the stationary-phase 
method. Numerical results illustrating the theory are presented. 



1 Introduction 

Oscillatory interlayer exchange coupling (IEC) has been found in a number 
of ferromagnetic/non-magnetic multilayer systems and is in some cases ac- 
companied by an oscillatory magnetoresistance. The physical origin of such 
oscillations is attributed to quantum interferences due to spin-dependent con- 
finement of the electrons in the spacer. The periods of the oscillations with 
respect to the spacer thickness can be correlated to the spacer Fermi surface, 
a relation frequently used in experimental studies. A number of models was 
proposed to explain this phenomenon and we refer the reader to excellent 
recent reviews on the subject 

The situation is much less satisfactory if the amplitudes and/or phases are 
concerned. They both depend sensitively on the details of the Fermi surface, 
and, from the experimental point of view, on the quality of the multilayers. 
Typically, samples include various amounts of disorder at interfaces as well 
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as in the bulk (e.g., surface roughness, intermixing, impurities, grain bound- 
aries, etc.) which can influence the amplitudes and the phases significantly. 
From the theoretical standpoint of view it is important to keep in mind that 
the IEC is an oscillatory phenomenon for which, strictly speaking, ampli- 
tudes and/or phases are defined only in the asymptotic limit. Experimental 
data, however, are usually only available for the first few oscillations which 
are sufficient to extract periods, but not amplitudes and phases, in particular 
for the so-called long-period oscillations. The presence of impurities not only 
complicates the theoretical studies but also can provide a valuable insight 
into the effects controlling the IEC. In particular, substitutional alloying can 
provide a valuable informations concerning the topology of alloy Fermi sur- 
faces. Alloying has also another, more subtle effect, namely it influences both 
amplitudes and phases and it can even introduce an extra damping of the 
oscillation amplitude (an exponential damping in addition to the usual 1 /N 2 
decay, where N is the spacer thickness) if ky-resolved electron states in the 
neighborhood of so-called callipers (extremal vectors of the Fermi surface) 
are influenced by disorder. Finally, we mention that a special case of alloying 
is intermixing of magnetic and spacer atoms at interfaces which can sig- 
nificantly influence coupling amplitudes and which occurs frequently during 
sample preparation in actual experiments. 

It is thus obvious that the study of the effect of alloying on the periods, 
amplitudes, and phases of the IEC is an important issue which, however, is 
not properly reflected in the available literature. Conventional bandstructure 
methods are of limited use for such studies although in particular cases, when 
combined with the virtual-crystal- type approximations (VCA), they may be 
justified, e.g., for VCr or CrMn alloy spacers studied recently [[!]. However, 
the complete neglect of alloy disorder makes a reliable determination of the 
coupling amplitudes or phases and, to some extent, even of the coupling 
periods, uncertain even in such favorable cases. 

In addition, reliable conclusions and verifications of experimental mea- 
surements can only be based on a parameter-free theory. In order to de- 
termine the IEC one typically estimates the energy difference between the 
ferromagnetic (F) and antiferromagnetic (AF) alignment of a system consist- 
ing of two magnetic slabs separated by a non-magnetic spacer. Using total 
energy differences (evaluated with the local density approximation to the 
density functional theory) represents an extremely difficult task as the tiny 
exchange energies have to be subtracted from the background of huge to- 
tal energies. Even if one employs very fast and accurate linear methods and 
computational tricks, the spacer thickness for which the calculated IEC val- 
ues are reliable, is limited to about 20 layers [pP]. On the other hand, for 
thin spacers this is the most accurate approach. One can alternatively employ 
asymptotic theories which are, strictly speaking, valid in the opposite regime, 
namely, for large spacer and magnetic slab thicknesses. The idea is to deter- 
mine reflection (transmission) coefficients for an isolated interface between 
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magnetic and spacer metals and the extremal vectors of the spacer Fermi 
surface. The former quantities then determine the coupling amplitudes and 
phases while the latter quantities their periods. In this case the calculations 
can be performed by using conventional bandstructure methods and, in addi- 
tion, they will provide a deep insight into the physical nature of the IEC j7j. 
Note, however, that neither of the above techniques can be extended to treat 
disorder nor can they be used to interpolate between two limits, namely, the 
case of thin spacers (preasymptotic region) and of thick spacers (asymptotic 
limit). For this a theory is needed which can bridge both the preasymptotic 
and the asymptotic region within a unified framework: IEC values for a large 
set of spacer thicknesses (say, for 1-100 atomic layers) can be analyzed in 
terms of discrete Fourier transformation in order to reliably determine not 
only periods, but also coupling amplitudes and phases. In addition, one can 
sample various subsets in order to analyze both the preasymptotic and the 
asymptotic regime as well as long-period oscillations. 

The basic idea is to determine the IEC directly by employing the so-called 
magnetic force theorem [|],[| for rotations in spin space rather then shifting 
atoms as in the conventional force theorem j[o|. We can thus use the same 
potentials for both the F and AF (or, in general, rotated) alignments of the 
magnetic slabs (the frozen-potential approximation) and consider only the 
single-particle (Kohn-Sham) energies. 

This allows a direct formulation of interlayer exchange coupling based on 
an application of the Lloyd formula |ll[] in order to evaluate the difference 
between the grand canonical potentials of the F and AF alignment. The first 
calculations of that type were performed by Dederichs's group in Jiilich fl^| . 
The method used in the present paper extends the above approach in three 
relevant aspects: (i) a reformulation within the framework of a surface Green 
function technique by which linear scaling of the numerical effort with re- 
spect to the number of layer s |T^ , P^| is achieved; (ii) a proof of the so-called 
vertex-cancellation theorem [|15| in order to study the influence of alloy dis- 
order on the properties of the IEC, and (iii) an efficient method for a fast 
and accurate evaluation of integrals involving the Fermi-Dirac distribution 
function in order to study effects of finite temperature [|l6|,[l7| . In the present 
paper we will review these particular techniques that were developed in the 
past few years and subsequently applied to a number of cases including al- 
loy disorder [ p^l9| , |20| , ^l| , p2[ . In addition, we have studied systematically the 
effect of non-magnetic cap-layers [^3[|4| on the periods, the amplitudes, and 
the phases of the oscillations of the IEC. 



2 Formalism 

In this Section we derive an expression for the IEC for in general non- 
collinearly aligned magnetic slabs embedded in a non-magnetic spacer. 
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2.1 Geometry of the system 

The system considered consists of a stack of layers, namely from the left 
to the right: (i) a semi-infinite (nonmagnetic) substrate, (ii) a left ferromag- 
netic slab of thickness M (in monolayers, MLs), (iii) a nonmagnetic spacer of 
thickness N, (iv) a right ferromagnetic slab of thickness M', and (v) a semi- 
infinite (nonmagnetic) substrate. The thickness of the ferromagnetic slabs 
may extend to infinity. Eventually, one of the semi-infinite substrates may 
be substituted by a finite nonmagnetic cap of thickness P interfacing semi- 
infinite vacuum. In general, the various parts of the system can consist of 
different metals, including disordered substitutional alloys. We assume that 
the spin orientation of the right magnetic slab is rotated by an angle 9 with 
respect to that of the left magnetic slab. In particular, the cases 9 = and 
9 = it correspond to the ferromagnetic and antiferromagnetic alignments of 
magnetic moments of two subsystems, respectively. 



2.2 Electronic structure of the system 

The electronic structure of the multilayer is described by means of the tight- 
binding linear- muffin tin orbital (TB-LMTO) method [p5f . In particular we 
employ the all-electron scalar-relativistic version as generalized to the case of 
random alloys, their surfaces and interfaces p6|]2l| . The key quantity of the 
formalism, the physical Green function G(z), is expressed via the auxiliary 
Green function g a (z) in the screened tight-binding LMTO representation a 
as 

G(z) = X a (z) + f i a (z)g a (z) f i a (z), (1) 
where 

g a (z) = (P a (z)-S a )- 1 . (2) 

Here S a is a matrix of screened structure constants S^ L R , L i, and P a (z) is a 
site-diagonal matrix of potential functions P^(z). The potential functions 
are diagonal with respect to the angular momentum index L = (£m) and the 
spin index a =\,\ while the structure constants are spin-independent. The 
potential functions can be expressed via the so-called potential parameters 
C, A, and 7 in the following manner 

where for matters of simplicity all indices are dropped. Similarly, the quan- 
tities A" and fi a in ([!]) can be expressed as 
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As only the screened representation will be used the superscript a is omitted 
in the following. 

A separate problem is the determination of potential functions P(z) for 
a given layered structure. Here we only mention that by employing the mag- 
netic force theorem we can use the same potential functions for the ferromag- 
netic and rotated (or, antiferromagnetic) alignments. For random systems 
treated within the so-called coherent potential approximation (CPA) the po- 
tential function P(z) is substituted by its coherent potential counterpart, 
V(z\ whereby the formal structure of the Green function (|^) remains the 
unchanged. The methods of determination of (coherent) potential functions 
for collinear alignments of magnetic moments in the present context can be 
found elsewhere p7| , p6[ . 



2.3 Definition of the IEC 



The exchange coupling energy £ x , evaluated in the framework of the magnetic 
force theorem, is defined as the difference of the grand canonical potential f2\ 
between the ferromagnetic (A = F) and antiferromagnetic (A = AF) align- 
ments of two subsystems, i.e. £ x = £Iaf — ^f- More generally, the quantity 
of the physical interest is the difference of the grand canonical potentials be- 
tween a rotated (6 ^ 0) and the ferromagnetic (8 = 0) alignment of the two 
magnetic slabs, namely, £ x (9) = 5Q(6) = Q{6) — 47(0). 

The grand canonical potential Q of a system is defined by 

/oo 
f(E,T,n)N(E)dE, (5) 

where N(E) is the integrated valence density of states, f(E, T, /i) is the Fermi- 
Dirac distribution function at the temperature T and the chemical potential 
/i of electrons. It should be noted that at zero temperature the chemical 
potential coincides with the Fermi energy Ep of the system. The integrated 
valence density of states is then given by 

N(E)=--1mf Tr G(E' + iO) dE' , (6) 

where Tr means the trace over lattice sites R, angular momentum indices 
L = (£m) and spin indices a. Using (HH), the following identities can be 
verified 

± \{z) = -A 2 (z) , ±P( z ) = ^(z). (7) 
az az 

Together with formula we find 

^ TrlnA(z) + Trln.g(z) = -Tr G(z) . (8) 
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The grandcanonical potential (||) is then expressed as 
1 f°° 

Q(T, fi) = Im / f(E, T, n) Tr In \{E + iO) dE 

T J-oo 

1 r 00 

Im / f(E,T,n)Tr]ng(E + iO)dE. (9) 

J-oo 

The formula in (JsJ) is the expression for the grandcanonical potential within 
the TB-LMTO method and for finite temperatures. 

The rotated magnetic configuration is characterized by the set of rotation 
angles — {9r} for all the lattice sites. In the reference (F) state all the 
angles Or — while in the rotated state 9r = 9 in the rotated magnetic layer 
and Or = for all other lattice sites. The quantities X(0,z) and g(0,z) for 
the rotated system are given by 

A(<9, z) = U(<9)A(0, z)U t (0) , g(0, z) = [U(0)P(O, z)^(0) - S]- 1 . (10) 

Here [U((9)]rr/ = 5rr>V '(Or) is the rotation matrix for spin 1/2 particles 
defined in terms of the single-site matrices U (Or) |29| 

uw=(_ c s (id 

where c = cos(6/2), s = sm(6/2),U(6) W (6)=V^ (6)V(9) = 1, and detU(0) = 
detU^(^) = 1. We note that in the rotated magnetic configuration P(0, z) = 
(J(0)P(O, z)U'(0) is generally a non-diagonal matrix with respect to the spin 
indices a, a'. 

The first term in (||) is independent of 9 because \(z) is site (and layer-) 
diagonal, it therefore does not contribute to the exchange energy £ x (9), i.e., 
it is sufficient to consider the second part only, 

i r°° 

f2(6, T,n) = Im / f(E, T, /i) Tr In g(9, E + iO) dE . (12) 

7f J —oo 

It should be noted that the above expression is valid only in the absence of 
spin-orbit coupling. 

The magnetic force theorem used here for the evaluation of the IEC was 
used also in related problems, e.g., for the evaluation of the exchange energies 
of two impurities embedded in a nonmagnetic host Q and then extended to 
the case of Heisenberg exchange parameters between two sites in a magnetic 
material ||. In the latter case the magnetic force theorem is valid only for the 
infinitesimal rotations while in the former case it is valid also for 9 = tt [BOl. 



2.4 Configurational averaging 

Keeping in mind applications to random systems, one is interested in the 
configurational average of the expression in ([l2]) , namely, 

(0)=--Im[ f(E,T, t i)(Trlng(E + iO))dE, (13) 
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where (. . .) denotes a configurational average. Difficulties here arise from the 
fact that the configurational average of the logarithm (In g(z)) can differ 
significantly from the logarithm of the configuration average \n(g(z)). The 
difference X = (In g) — In (g) , the so-called vertex correction, is difficult to 
calculate and usually cannot be neglected. Fortunately, this problem can be 
circumvented by using the vertex cancellation theorem [ p^[ , which states that 
the contributions from the vertex correction for the F and AF configurations 
cancel each other exactly, namely Tr Xaf~ Tr X^ = 0, such that to first order 
with respect to the angle between the magnetizations in the two ferromagnetic 
layers vertex corrections can be omitted. In other words, the evaluation of 
(|13|) simplifies to 



(/?) = --Im f f(E, T, (j.) Tr In (g(E + iO) > dE , 

71 J-oo 



1 Im / f(z, T, n) Tr In (g(z) ) dz . (14) 



7T 



We have also substituted the energy integral by integration over a contour 
in the complex energy plane z. The possibility to neglect vertex corrections 
can conveniently be used in calculations of the interlayer exchange coupling 
as explicit numerical calculations have shown that it remains valid to a good 



accuracy even for an angle as large as 7r 15 . In this respect it is very sim- 
ilar to the force theorem pp| . It is important to note that such an exten- 
sion is only applicable to the evaluation of exchange energies of magnetic 
systems interacting via a non-magnetic host. An evaluation of exchange en- 
ergies in ferromagnetic systems such as parameters of a classical Heisenberg 
model, was claimed to be limited to infinitesimal rotations only ||. The use 
of the vertex-cancellation theorem allows to reduce the computational time 
in first-principles calculations by almost two orders of magnitude, so that 
the computational effort for disordered systems is comparable to that for a 
pure system jL5). We refer the reader to Appendix A for more details con- 
cerning the derivation and applicability of the vertex-cancellation theorem. 
The last remark concerns the fact that the expression for the change in the 
grandcanonical potential within the magnetic force theorem also includes the 
classical magnetostatic dipole-dipole interaction energy (DDIE). The DDIE 
decays with a spacer thickness much faster than the IEC and its contribution 
can be thus neglected for thicker spacer anyhow. In addition, first-principles 
fully-relativistic calculations of the IEC |3^] have demonstrated that this term 
has a negligible influence even for a rather thin spacer amounting just to a 
few layers. Consequently, the DDIE term will be neglected in the following. 



2.5 Lloyd formula 

We need to evaluate the difference of configurationally averaged grandcanon- 
ical potentials in the rotated and FM configurations. This can be done conve- 
niently with the help of the well-known Lloyd formula p] applied to layered 
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systems. We formally split the system into two non-interacting fragments, 
namely a left fragment C, which consists of the left substrate and the left 
magnetic slab, and a right fragment 1Z, which comprises the rest of the sys- 
tem, i.e., the spacer, the right magnetic slab, and the right substrate (or, 
eventually, the cap layer interfacing the vacuum). Fragments are described 
by the unperturbed Green function (go(z)}. In the next step we couple two 
fragments together with help of a localized potential V which is simply the 
intcrlayer screened structure constant. This procedure has a number of ad- 
vantages as compared to a conventional way of embedding two finite mag- 
netic layers into the infinite (bulk) host spacer Jl2|: (i) the perturbation V 
is independent of the thicknesses of magnetic layers; (ii) complicated sam- 
ple geometries can be treated, including semi-infinite magnetic layers; and 
(iii) a powerful and efficient method exists for the evaluation of the Green 
function of fragments, namely the surface Green function technique in the 
principal-layer formulation |p6| , p7| . 

Keeping in mind the vertex cancellation theorem, one gets for a difference 
in the configurationally averaged grandcanonical potential (|l4l), the expres- 
sion 

(Sf2) = --lm[ f(z,T,v)Tr\n(l-V{g (z)))dz, (15) 

7T JC 

where (go(z)} is the configurationally averaged Green function of the de- 
coupled non-interacting fragments C and 1Z defined above. For the sake of 
simplicity, we will denote from here on the configurationally averaged quanti- 
ties by an overbar, e.g., (go{z)} = go(z). The concept of principal layers (PL) 
[ p3| as used within the TB-LMTO method leads to a block tridiagonal form 
of the structure constants and of the inverse Green function. If we apply this 
tridiagonality to (|l5|), we get for V and (go(z)} the following expressions by 
using a supermatrix notation with respect to nearest-neighbor PLs resolved 
in the wave- vector k|| , 



,(16) 



where 5io(k||) = [5oi(k|i)l . Combining ( p"5| ) and ( |l6| ) one gets 

5Trlng(z) = -^^trln [l-r £ (k||,z)^(k||,z)] , (17) 
II k y 

f iC (k||,«) = 5io(ky)gc(k|| ) z)5oi(k||). 

Here the quantity /!c(k|| , z ) nas the meaning of an effective embedding poten- 
tial, and the quantities Qc and Qn are the configurationally averaged surface 
Green functions (SGF) Q of the magnetic subsystems C and 7Z, respec- 
tively. By definition, the surface Green function Q$ (S = C, TV) is the top PL 
projection of the Green function of the corresponding semi-infinite system S. 
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Its determination in the case of random systems was extensively discussed in 
the literature, see ]3^ , |35| , ^6p^ | . The summation in ([n]) extends over the sur- 
face Brillouin zone (SBZ) corresponding to the underlying two-dimensional 
translational symmetry j3?j], and iV|| is the number of sites in a layer. 

2.6 The IEC for a general angle 6 

Let us now turn to the evaluation of the energy difference between arbitrary 
alignments. Consider the following quantity, 

tr In Z = tr In (1 - A B) - tr In (1 - A B ) , (18) 

where the matrices Ao and Bo are related to the ferromagnetic alignment and 
thus are diagonal in spin space 

'kl A „ fBl 

The particular form of the subblocks Aq and Bq (a =f, J.) is given by 

AJ = S 10 (k||)52(k||,z)5oi(k||), B^^(k,|,z). (20) 

The matrix B refers to an alignment in which the orientations of the 
magnetization in two magnetic slabs are rotated uniformly by a relative angle 
0, 

B = U(0)B o Ut(0), (21) 



A ° = ' o°aJ' Bo = obV' (19) 



where U(0) is the rotation matrix (|l l|). The quantity 1 — Ao B in ( |18[) can 
therefore be written as 

1 - A B = (\J(6) - A o U(0)B o ) U* (9) , (22) 
where, as follows from (B~9|) and (0), 



U(0) - Ao U(0) Bo = ( C( ;\ A i B ^ 'I !! t? • (23) 



vl Bj) s (1 - Aj B* 
-a (1 - A l Q B ) c (1 - A b£) 

Using now the identity trln X = lndet X, which is valid for any non-singular 
matrix X, and the identity 

dct ( C D ) = dctA • detD ' det ( 1 _ A 1 B D 1 C ) ' ^ 

which in turn is valid, if the matrices A and D are non-singular, it is straight- 
forward to prove that 

/ 1 - cos(0) \ , . 

trlnZ = tr L ln 1 — — M , (25) 
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where 

M = 1 - (1 - A] Blr 1 (1 - A B l Q ) (1 - Al Bl)- 1 (1 - A l B ) . (26) 



It should be noted that in (|18|) tr denotes the trace over angular momenta 
and spin, while in (|24|) tr^ denotes the trace over orbital momenta only. The 
final expression for £ x (9) is thus given by 

tr L In (l - L^g(g) M (k|| ,z))dz, (27) 

in which the energy integral is expressed in terms of a contour integral which 
will be discussed in detail later. 

It is interesting to note that the expression (26) for M(kii, z) can be rear- 
ranged in the following form [jl9| 

M = - fl - Sxo G\ Soi Ql) ~* S 10 (§1 - G l c 



1 - Soi G^k Sio G l c) Soi [Gk - G l n ) ■ (28) 

It explicitly factorizes the 'spin-asymmetry' of the problem and it is directly 
related to RKKY-like theories Q. This result [|l9) is formally equivalent to 
the results of the spin current approach fl39| as formulated within a Green 
function formalism based on an empirical single orbital tight-binding model 
[ [l0| . A matrix version developed in the framework of a semiempirical tight- 
binding model has appeared recently |^] . 

For completeness we also give the result for the common case of the an- 
tiferromagnetic alignment (6 = tt): 



S x = S x (n) = ^-J2 lm I /(«,T,/i)tr £ lnA^(k||,«)cUr, (29) 

where M. is a product of four terms, 

M = (1 - Al Bl)- 1 (1 - Al B l Q ) (1 - Al Bl)- 1 (1 - A l Q Bj) . (30) 



2.7 The torque and infinitesimal rotations 

The differential change in the grand canonical potential 5f2(8) with respect 
to a differential relative angle 9, —dS(l(9)/d9, is usually called the torque. 
The torque can easily be obtained by differentiating ( |27| ) with respect to the 
angle 9. By definition one gets therefore 

T{6) = -^fi or £ X (6)=-J T(6')d9', (31) 
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whereby T{9) follows immediately from (]2 



tr L 



Ulk .:) | I -i[l-oos(0)]M(k||,*; 



-l 



dz. (32) 



By formally expanding the logarithm in ( |27| ) in powers of 1 — cos(#), one can 
cast the expression for £ x (6) into the form 

£ x {9) =B x [l- cos(0)] + \b 2 [1- cos(0)] 2 + . . . , (33) 

where B\ and B2 are the so-called bilinear and the (intrinsic) biquadratic 
exchange coupling coefficients, respectively, 

Bl = YJT E Im / T ' tr ^ M ( k H > z ) dz > ( 34 ) 

S 2 = £ Im / /(*, T, M ) trz [M(k,| , z)] 2 dz . 

^ 11 k n C 

It may be, however, more convenient to fit the exact expression (E7|) into the 
form d33| ) by employing calculated values for 8 = tt/2 and 6 = n Q. We 
obtain 

g a (7r) + 2g a (7r/2) £ x (n)-2£ x (n/2) 
B\ = 5 ' B2 = 2 ' ^ ' 

Of particular interest is the expansion of £ x (0) for a small 9, i.e., when 
1 — cos(#) is a small parameter (the method of infinitesimal rotations (MIR)). 
This approach becomes particularly relevant in the case when the spacer is 
a magnetic metal or for complicated geometries, e.g., for so-called periodic 
multilayers. 



2.8 The IEC as interface-interface interaction 

We will now discuss briefly an alternative approach of a direct evaluation of 
the IEC as a difference in the interface-interface interaction energies rather 
then its indirect determination in terms of the energy of a single interface ( |l3| - 
[l6| ). We decouple the system into three fragments, a left, central, and right 
fragment, £, C, and 1Z, respectively. The left and the right fragment are 
formed by corresponding substrates with magnetic slabs whereby the central 
slab comprises the spacer. Both approaches are physically equivalent because 
it is irrelevant how the system is divided into an unperturbed part and a 
perturbation. Note, however, that the interface-interface formulation is more 
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general as it could be used for a determination of interaction energies of two 
generally different interfaces. 

The derivation proceeds in two steps and employs partitioning technique 
with respect to the trace of the logarithm of the Green function. First, the 
subsystems C and 1Z are downfolded which leads to an effective problem 
of two localized perturbations in the subsystem C. The second step, a two- 
potential formula applied to the fragment C separates directly the interface- 
interface contribution. The result has formally the same structure as the 
previous one ( [l7| , f26"l ) , but the subblocks Aq and Bg (cr =|, J.) are now of the 
following form 

A-o = <7ivi(k||, z) rf (k|| , z) <hiv(k|| , z) , Bq =r^(ky,z) . (36) 

The t- matrices fj (i — 1,N) corresponding to "multiple scattering" at indi- 
vidual interfaces C/C, (i — 1) and C/TZ, (i — N) are expressed as 

7f(ky,z) = /7(k||,*) [l-gu^ h z)rr^ h z)Y 1 , (37) 

where the effective embedding potentials /^(kiuz) of the left and right in- 
terfaces (i = 1,N), respectively, are defined as 

rf (k|, , z) = S l0 (k,| ) Ql (k,| , z) S 01 (k,| ) , (38) 
r^(k|| J «) = S 01 (k||)^(k|| J z)Sio(k||). 

Here, Qg (S = £, TZ) are the configurationally averaged SGFs of the left and 
the right semi-infinite regions, respectively. Details of the derivation can be 
found in Appendices B and C. The coupling between the two magnetic sub- 
systems is due to the layer off-diagonal projections </ijv(kii , z) and gjvi(kii , z) 
of the Green function (GF) of the finite spacer consisting of N layers. The 
oscillatory behavior of interlayer coupling is then governed by the oscillatory 
behavior of these quasi one-dimensional spacer Green functions, a formula- 
tion which is very much in the spirit of a simplified RKKY approach JlJ. 
An efficient method of evaluation of the corner-blocks of the Green function, 
<7ij(ki|, z), — 1, N), is described in Appendix D [^2|j36| . 



2.9 Relation to the KKR method 

We shall discuss now the relation of the present approach ( ^9| , ^o| , ^6|) to the 
method employed in jl2| and based on the Korringa-Kohn-Rostoker (KKR) 
Green function technique. Let us note first the deep internal connection be- 
tween the KKR and the TB-LMTO-GF approach (see §1,0 for more de- 



tails). The model in (12) consists of an infinite ideal non-magnetic spacer as 
a reference system and of two magnetic slabs representing localized pertur- 
bations. For simplicity we start from the case of two magnetic monolayers in 
an infinite spacer. The result 

A% = g b m (z)(k h z)q(k h z)g b 1N (k h z), B° = t%(k h z) (39) 
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is formally the same with the exception that the r-matrices entering (|3^) are 
now substituted by the single-site t-matrices U which describe the scattering 
of electrons from two magnetic monolayers at i = 1, N embedded in an infinite 
non-random bulk spacer and separated by N — 2 spacer layers: 

%. i H ,z) = APZ. i {z) [l + g b (k h z)APZJz)]- 1 . (40) 

The strength of the scattering potential, APf -(z), is given by the difference 
of the potential functions for the magnetic monolayer Py^z) and for the 
non-magnetic spacer P(z), while g b (k\\,z) is the layer diagonal block of the 
GF of the bulk spacer. The layer off-diagonal blocks of the bulk spacer GF, 
G\ N {z) and G h N1 {z), are given by 

g\ N m,z) = [g s m,z)s 01 {^] N - 1 g b (k h z), (4i) 

and similarly for Q b Nl (z). Here, !? s (k|| , z) is the corresponding SGF of an ideal 
semi-infinite non-magnetic bulk spacer |53| . It should be noted that also the 
layer-resolved bulk Green function S b (k|| , z) can be expressed in terms of the 
SGFs (see, e.g., Since ( [flf ) is exact, there is no need to perform an 

additional fc^-integration p^ |. It is easy to show that the result is formally 
identical to the case of two impurities in a simple tight-binding linear chain 
model with nearest neighbor hopping. 

A generalization to the case of magnetic slabs containing a finite number 
M of magnetic layers is formally straightforward The t-matrices t^.Az) 
are then supermatrices with respect to angular momentum and layer indices 
and the numerical effort to evaluate (^o| ) increases with the third power of M 
as contrasted with the results of the present approach (|l^,^6|) which depend 
only linearly on M. 



2.10 Influence of external periodicity 

Until now it was assumed implicitly that we have a simple "parent" lattice 
[ [37| . The periods of the coupling oscillations are closely related to the Fermi 
surface geometry jjJU of the bulk spacer. A different translational symme- 
try (complex lattices) or stacking sequence within layers will thus tend for 
sufficiently thick spacers to a different kind of bulk periodicity and hence to 
new periods. For example, an alternating stacking of fcc(001)-layers Cu and 
ordered c(2 x 2)-CuAu layers tends to an ordered fcc-Cu3Au alloy with a 
Fermi surface topology different from that of fcc-Cu spacer. For a discussion 
of " superlattice" formation in magnetic multilayers see also |38| . We will dis- 
cuss in the following in some detail two possibilities, namely superstructures 
in the spacer and in the magnetic slabs. 

We start with the former case by assuming the same geometry as discussed 



in Sec. 2.1 but now the spacer slab consists of two non-magnetic metals A and 
B with respective thickness tla and Ub periodically alternating. Typically, 
the spacer layer starts with the layer A{B) and ends with the layer B(A), 
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but the termination of the spacer slab with the same layers is also possible 
(and interesting p2|]). The particular case of ua = Ub = 1 corresponding to 
an (OOl)-stacking of an ordered fcc-CuAu alloy was already treated on a first- 



principles level 1 22 . The more general case, (nA,ns > 1), which corresponds 
to artificially grown superstructures, was treated only within a simple one- 
band model p4| . In both cases, new periodicities (in comparison with the 
spacers consisting from pure A or B metals) arise with an increasing number 
of repetitions. Alternatively, one can consider a superstructure within a given 
spacer layer, or combination of both, e.g., the above mentioned example of 
the ordered fec-CuaAu alloy spacer. The similar situation can be encountered 
also in the magnetic slabs. In particular the case of a c(2 x 2)-CoFe periodicity 
within the magnetic layers separated by a fcc-Cu(OOl) spacer |^(| leads to 
the rather surprising appearance of new periods. These new periods can be 
now correlated to critical points of the spacer Fermi surface folded down to 
the Brillouin zone corresponding to a c(2x2)-superlattice |Q. A correlated 
gradual appearance of new periods and the order in statistically disordered 
layers is a clear indication of their relation to a different bulk periodicity 

Hi- 

A special case of alternating layers of A and B metals is when one of 
metals is magnetic and the other is nonmagnetic, all of which sandwiched 
between two substrates. This is the case of a periodic multilayer. 

The generalization of the present formalism to above discussed cases is 
rather straightforward. In the case of a superlattice within a layer it is just 
sufficient to substitute matrices appearing in ( p7| , p9| ) by the corresponding su- 
permatrices, e.g., by (2x2)-supermatrices in the case of a c(2x2)-superlattice. 
The key quantity, the surface Green functions K ( pp| ) , can be easily evalu- 
ated also in this case (see for details ^f|). The generalization of the formalism 
to the case of alternating layers from A and B metals is as well simple be- 
cause the surface Green function is constructed in an epitaxial manner, i.e., 
layer by layer, and it is therefore immaterial if the stacking of layers consists 
of the same or a different material. In the limit of a periodic multilayer we 
should just keep in mind that a proper repeating unit consists now from four 
layers, namely S — M — S — M, where the symbols S and M refer to the 
spacer and magnetic layers, respectively. This is necessary to calculate the F 
and AF configurations needed for the evaluation of the IEC. We note that 
the present formalism allows to evaluate efficiently and reliably the IEC for 
thick spacers (one hundred layers and more) which is important for realistic 
studies of so-called superlattice spacers and of periodic multilayers. 



2.11 Temperature-dependence of the IEC 

We conclude this Section by reviewing a recently developed technique for an 
efficient evaluation of the temperature dependence of the IEC (l6) . The main 
cause for the temperature dependence of the IEC is connected with thermal 
excitations of electron-hole pairs across the Fermi level as described by the 
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Fermi-Dirac function. It turns out that other mechanisms (as for example 
electron-phonon and electron-magnon interactions) are less important. We 



rewrite (29) into the following form 

£ x (T) = lmI(T), J(T)= f f(z,T,p)¥{z)dz, (42) 
where 



c 



n 11 k n 

with the energy integration performed over a contour C along the real axis 
and closed by a large semicircle in the upper half of the complex energy plane. 

The integral in ( J42| ) can be recast into a more suitable form using the 
analytic properties of ^(z), namely, (i) ^(z) is holomorphic in the upper half 
of the complex halfplane, and (ii) z^P(z) — > for z — > oo, Imz > 0. Let us 
define a new function <P(y) = —i <P(Ef + iy) of a real variable y, y > 0. Then 
at T = K, 

r + CO 

1(0) = / $(y)dy, (44) 
Jo 

while at T > K, 

OO 

I(T) = 2nk B Tj2$(yk), (45) 
k=l 

where kg is the Boltzmann constant and the yk are Matsubara energies, 
y k = Trk B T(2k - 1). In the limit T — * 0, I(T) -> 1(0) continuously. 

We have verified that the function <I>(y) can be represented accurately as 
a sum of a few complex exponentials of the form 

M 

${y)= y £ J A j exp(p j y), (46) 
j=i 

where the Aj are complex amplitudes and the pj are complex wave numbers. 
An efficient method of finding the parameters Aj and pj is described elsewhere 
[ fl6| . The evaluation of I(T) is then straightforward: 

M . 

I(T) = -2vk B T £ i , (47) 

•f^ exp (7tkb J j?j) — exp (— 7tkbJ pj) 

which for T = K gives 

J(0) = -E£. ( 48 ) 
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The efficiency of the present approach allows to perform calculations with a 
large number of k|| -points in the irreducible part of the surface Brillouin zone 
(ISBZ) in order to obtain accurate and reliable results. Note also that such 
calculations have to be done only once and then the evaluation of the IEC 
for any reasonable temperature is an easy task. 

The effect of finite temperatures on the IEC can be evaluated also ana- 
lytically. The analytical approach assumes the limit of large spacer thickness, 
for which all the oscillatory contributions to the energy integral cancel out 
with exception of those at the Fermi energy. The energy integral is then eval- 
uated by a standard saddle-point method Q] . The general functional form of 
the temperature-dependence of the interlayer exchange coupling £ x (T) in the 
limit of a single period is then given by 

cNT 

S x (T) = S x (0)t(N,T), t(N,T)= . (49) 

smh(cA< 1 ) 

Here, N denotes the spacer thickness in monolayers, and c is a constant which 
depends on the spacer Fermi surface. The term £ x (0) exhibits a standard 
iV~ 2 -dependence [|l], while the scaling factor t(N, T) depends on the product 
N and T. In the preasymptotic regime (small spacer thickness) the functional 
form of t(N, T) differs from that of (^9|), particularly in the case of a complete, 
but relatively weak confinement due to the rapid variation of the phase of the 
integrand which enters the expression for the IEC |Q . The present numerical 
technique is free of the above discussed limitations and can be used to check 
conclusions of model theories. 



3 Numerical results and discussion 
3.1 Details of calculations 

Special care has to be devoted to the energy and the Brillouin zone integra- 
tions. For a finite temperature we determine the parameters of the complex 
exponentials in ([l6]) through an evaluation of ^{y) at 40 Matsubara ener- 
gies corresponding to T = 25 K. We have verified that the results depend 
weakly on the actual value of the parameter T. For T = K we have tested 
two energy contours C, namely a semicircle between the bottom of the band 
(E m i n ) and Ep, or, alternatively, a line contour Ep + ie, e G (0, oo), using 
a Gaussian quadrature. The results were very similar in both cases. Using a 
line contour avoids possible problems connected with the phase of a complex 
logarithm. Typically a total of 10-15 energy points was used. A large number 
of ky-points in the ISBZ is needed only for energy points close to the real 
axis, whereby generally a greater number is needed for lower temperatures 
and thicker spacers. The number of k|| -points can significantly be reduced 
for energies well off the real axis. In particular, for the first energy point on 
the contour close to the Fermi energy we typically use 5000-10000 k||-points 
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in the ISBZ, while for the next 3-4 energy points the number of k|| -points 
is reduced by a factor two for each other point, and about 50-100 kii-points 
are taken for all remaining energy points on the contour. The thickness of 
the spacer, for which well converged results are obtained, is about 100 spacer 
layers. 



3.2 Analysis of the results 

The calculated results, namely £ x (0, N), where N specifies the spacer thick- 
ness, can be analyzed in terms of a discrete Fourier transformation 



N max 

F(6,q) = ~ V N 2 £ X (6,N) exp(iqN) , (50) 

y N=N min 

where p = N max — N min + 1 is the number of values used in the Fourier anal- 
ysis, and N min is chosen in order to eliminate the effect of very thin spacers, 
or, to analyze intentionally either the preasymptotic or the asymptotic re- 
gion. Typically p is about 40. The background oscillations thus occurring [jl4| 
are due to the discreteness of the Fourier transformation. The background 
oscillations can be smoothened using the procedure described in Q , namely 
by multiplying N 2 £ x (9, N) by C snv(nN / p) / (ttN / p) , where C is a normaliza- 
tion factor. The periods of oscillations A a (in monolayers) arc then identified 
with the positions q a of pronounced peaks of ^(g^)! as A a — 2ir/q ai the 
amplitudes of oscillations A a are estimated from A a = (2/p)\F(q a )\, and 
their phases from <f> a — ir/2 — AigF(q Q ), (a = 1, 2, . . .). This analysis can be 
extended to more complicated cases, namely when the IEC is a function of 
two variables, e.g., as a function of the spacer and cap thicknesses N and P, 
respectively. A two-dimensional discrete Fourier transformation 

N 2 P 2 

F2(0,q N ,q P )= £ ^ (N + P) 2 £2(8, N, P) S qNN+qpP ^ (51) 
N—Ni P=Pi 

is a suitable tool to analyze the quantity £2(0, N,P), where the prefactor 
(N + P) 2 is consistent with the asymptotic behavior [|23|j47t| for large spacer 
and cap thickness. Strictly speaking, this is quite an obvious choice for the 
case when the spacer and cap are formed by the same material, but it can be 
used also when the spacer and the cap correspond to different materials (for 
more details, see p3fl ). In ( |5l|) we have introduced the quantity 

£ 2 (9,N,P) = £ X (9,N,P)-£ (9,N), £ Q (0,N)= lim £ x (0,N,P) (52) 

P— >oo 

in order to remove a trivial peak in the absolute value of F2{0,qN 1 qp) at 
Qn = Qp = 0. A similar two-dimensional discrete Fourier transformation 
is also useful in the study of the IEC with respect to the thicknesses of the 
spacer and the magnetic slabs. We note that if one of variables, e.g., the spacer 
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thickness N is fixed, it is possible to analyze the calculated IEC values again 



with the help of (50) 



An alternative of calculating the Fourier transform (50) consists in sub- 
dividing the k||-integral in ( p7| ) into areas around the critical k||-vcctors (cal- 
lipers) related to the different oscillation periods |3(^9|. In the asymptotic 
limit each subarea gives then rise to a single oscillation period, while in the 
preasymptotic regime the resulting division into different periods is only qual- 
itatively valid. In a sense this method bridges the present method of discrete 
Fourier transformations and the purely asymptotic trea tme nt of calculating 



only the behavior of the critical k|| -vectors (see Section p.3|) 



3.3 Asymptotic expansion 



Model studies [lj^] indicate that in the asymptotic region, i.e., for large spacer 
and magnetic layer thickness, and for a random sample, the general form of 
the spacer-thickness dependence of the IEC is given by 

7 

£ x = Im ^2 j£ exp(iQ a N) . (53) 

a 

Here the sum runs over all possible periods a, the quantities Z a and Q a are 
the complex amplitudes and complex stationary points (callipers), respec- 
tively, defined in the following manner 

Z a = A a exp(i^ Q ) , Q a =q a + i^a- (54) 

The quantities A a and <P a are the amplitudes and phases of coupling os- 
cillations, p a = 2tt / 'q a their periods, and the quantity A Q characterizes the 
damping of oscillations due to the effect of alloying in the sample determined 
at the Fermi energy. In the limit of non-random samples, A Q = 0. 

The parameters in (p3|) can be extracted from a detailed knowledge of the 
spacer Fermi surface J7[. We briefly sketch a numerical way of determining 
of the parameters of this asymptotic expansion which requires the knowledge 
of the integrand of ( |2S| ) for a set of k|| -points in the neighborhood of the 

stationary points k|" . 

The expression ( p9| ) for IEC at T — K can be rewritten as 

£ X = ^-Im V F(k,|), r(k|,) = i f f(z,0)tv L ]nM(k b z)dz. (55) 

N \\ t; n Jc 

The integration with respect to the energy variable is performed numerically. 
The function 5^(k||) for large N decreases as 0(1/ N) and behaves like 

F(k||) = ^expW(k||)), (56) 

where the pre-exponential factor <?(kii) is a smooth function of kn and the 
phase 4>(k\\ ) has one, or more stationary points in the SBZ that correspond 
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to callipers of the spacer Fermi surface such that Vk,, 0(kii) = 0. The integral 



over the SBZ in (£>5|) can be evaluated using the stationary-phase method. 
The contribution of a stationary point kj"^ = (kx , ky ) is found in the 
following way: the integration limits are extended to infinity, and the phase 
function </>(ky) is approximated by a quadratic function of k|| = (k xi k y ) in 
the vicinity of the stationary point, 



(a)x 



0(k,|) = 0(k 



E 



E 



Qij(ki k\ )(fc 



(57) 



The expansion coefficients Qy, Pi, and 0(kj| Q ^) are determined by a least- 
square fit to values of <f>(\t.\\ ) calculated in the vicinity of ky . This procedure 
allows to eliminate numerical inaccuracies with respect to both the values of 
Qij and the position of the stationary point k n , and it is applicable even 
for disordered surfaces. By inserting ( |56| ) and (57) into ( j55|) we find 



1 



7riVVsBZ 



Im { g(k^) x 



f 5(ki Q) ) 
Im \ t : = exp 

SBZ 



N 2 V<- 



V-detlQI 



zA^(kf, Q) ) 



(58) 



where the two-dimensional integration region D extends to infinity, and Vsbz 
denotes the volume of the SBZ. The second line in (|58|) is obtained by di- 
agonalizing the quadratic form in the exponent ([57]) and by evaluating the 
resulting one-dimensional Gaussian-like integrals. The identification of the 
parameters is now straightforward, namely 



fl(k[, a) ) 



^sbz ^/-dct|Q| ' 



Qa 



0(kf, tt) ). 



(59) 



3.4 Free-electron limit 

The numerical efficiency of the present formalism offers an interesting possi- 
bility of testing model theories p]]. The simplest of such models is the free- 
electron model, because of a spherical Fermi surface with a single critical 
vector at k|| = and a trivial correspondence between the value of the os- 
cillation period and the band-filling. The free-electron model can be easily 
simulated by the present formalism by replacing the true metallic potentials 
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by flat potentials (the empty-sphere model). In this case the potential func- 
tions (||) are analytical functions of the lattice constant. For a suitable choice 
of the lattice constant and the position of the Fermi energy it is irrelevant 
what lattice and layer stacking is used, e.g., the fcc(001)-stack is the simplest 
choice. On the other hand such a model is free of the limitations usually 
adopted jf], e.g., the assumption of large spacer and magnetic slabs thick- 
nesses, or the approximate evaluation of the energy integral for the case of 
finite temperatures. 



3.5 Numerical illustrations 

In Fig. 1 N 2 £ X (N) is displayed as a function of the spacer thickness N for two 
semi-infinite Co(001) subsystems sandwiching an fcc-Cu spacer. The corre- 
sponding discrete Fourier transformation in Fig. 2 shows a pronounced short- 
period oscillations of 2.53 monolayers (MLs) while the long-period oscillations 
are suppressed in this geometry ]l3| , p^|j4^ ] . The results are insensitive to the 
choice of the lower and upper index in the summation in (^0|) provided the 
preasymptotic region is excluded |H. 



20 n 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 r 




5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 
N (monolayer) 

Fig. 1. Exchange coupling N 2 £ X (N) at T = K as a function of the spacer thickness 
N for two semi-infinite fee Co(001) subsystems sandwiching a Cu spacer. Diamonds 
refer to the calculated values, the full line (back Fourier transform) serves as a guide 
to the eye 



For a large enough N the IEC can be approximated by the asymptotic 
form in (53). The amplitude, phase, and the wave- vector entering this expres- 
sion can be determined from the calculated £ x (N) in the manner as described 
in Sec. |3.2| and the asymptotic result (^3|) was compared with the calculated 
results for a large set of systems including both ideal and alloyed semi-infinite 
fcc(001) magnetic subsystems sandwiching a Cu-spacer: overall good agree- 
ment was found An example of the complex amplitude for this case 
is presented in Fig. 3 illustrating the insensitivity of the phase to elements 
which form the magnetic layers. It is seen that phases corresponding to Co, 



Interlayer Coupling 21 




Fig. 2. Absolute value of the discrete 
Fourier transformation of N 2 £ X (N) for 
a finite set of spacer layers (iV=20-80) 
corresponding to two semi-infinite fee 



12 3 Co(001) subsystems sandwiching a Cu 

q-vector spacer. The temperature is T = K 



Fe5oNi5o, and Fei/3Ni 1 / 3 Coi/3 which have the same average electron numbers 
N e i=9 are nearly the same Mj. 




Fig. 3. Complex amplitude Z 1 ^ 2 = 
A 1/2 e l4>/2 , where A and <P are the os- 
cillation amplitude and phase, respec- 
tively, for a semi-infinite fcc(OOl) sub- 
systems formed by Fe, Co, Ni, their bi- 
nary alloys (bullets), and the ternary 
alloy Fei/3Coi/3 Nii/3 (square) sand- 
wiching a Cu spacer. The units are 
(mRy) 1//2 . The dotted, dashed, and 
full lines connect various alloys and 
serve as a guide to the eye. The rays 
starting at the origin show approxi- 
mately the phase corresponding to the 
indicated average number of valence 
electrons 



The IEC depends on the temperature T via a factor x/sinh(x), x — cNT, 
where T is the temperature and N the spacer thickness. This remarkable 
result of model theories |jj was verified by calculations such as illustrated in 
Fig. 4. 
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1 1 1 1 1 1 

5 10 15 20 25 

NT/1000 (ML*K) 



Fig. 4. £ X (N,T)/£ X (N,T = 0) plotted 
as a function of f = AT for a trilayer 
consisting of semi-infinite fee Co(001)- 
slabs sandwiching a Cu spacer. The 
thick line refers to x/sinh(x), x = c NT 
with c = 0.000195 obtained by a least 
square-fit to the computed data 



The IEC depends in an oscillatory manner not only on the spacer thickness TV 
but as well on the thickness P of a covering cap. The oscillations are around 
a biased value which corresponds to coupling for a given spacer thickness 
assuming a semi- infinite cap. This phenomenon is illustrated in Fig. 5 in 
terms of discrete Fourier transformations with respect to the spacer and the 




0.0 



Fig. 5. Absolute values of the discrete two-dimensional Fourier transformation of 
(N + P) 2 £2 (AT, P) with respect to the spacer and the cap thickness in the case of 
two magnetic slabs each five monolayers thick with a Cu-substrate, a Cu-spacer, 
and a Cu-cap. For a definition of £i{N, P) see the text 



cap thickness (see Sec. [T^) for a sample consisting of a semi-infinite fec- 
Cu(001) substrate, left and right magnetic layers each five MLs thick, a spacer 
with varying thickness N, and a Cu-cap of varying thickness P. Fig. 5 shows: 
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(i) long-period oscillations (missing in Fig. 2) in addition to the short-period 
ones, and (ii) oscillations with respect to the cap thickness which are exactly 
the same as for the spacer because both are controlled by the same Fermi 
surface, namely that of fcc-Cu. The more complicated case of different spacer 
and cap materials is discussed in 23 24] . 




Fig. 6. Absolute values of the discrete 
Fourier transformation of N 2 £ x (N) 
for two semi-infinite fee Co5oFeso(001) 
subsystems sandwiching a Cu spacer 
with different kinds of chemical order 
in magnetic layers: (a) 3 = 1 (com- 
plete c(2 x 2)-order, full line), (b) S = 
0.8 (dashed line), (c) S = 0.5 (dashed- 
dotted line), and (d) S = 0.0 (disor- 
dered case, dotted line). The temper- 
ature is T = K 



Ordering in the spacer [|22| or in the magnetic layers [|20| can induce new 
periods due to the formation of two-dimensional sublattices. The situation is 
particularly interesting for a c(2 x 2)-ordering in magnetic layers sandwiching 
an ideal Cu-spacer p(J. As illustrated in Fig. 6 for full ordering two new 
periods with complementary periods and phases are formed in addition to 
a conventional short-period due to a fcc-Cu spacer pcj] . These new periods 
vanish in the completely disordered case. 

Finally, the effect of disorder in the spacer |lj| is illustrated in Fig. 7. 
Alloying of Cu with Ni decreases the number of average valence electrons and 
leads to a contraction of the alloy Fermi surface, and in turn to a reduction 
of the coupling oscillations. The opposite behavior has to be expected for 
alloying of Cu with Zn, whereas only a small concentration dependence of 
the periods for the CuAu case is seen. The amplitudes of the oscillations are 
generally reduced by alloying, and in the case of CuZn spacer they are even 
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5 10 15 20 25 30 35 40 45 



N (monolayer) 

Fig. 7. Exchange coupling N 2 £ X (N) at T = K as a function of the spacer thick- 
ness N for two semi-infinite fee Co(001) subsystems sandwiching a spacer of (from 
bottom to top) Cu, CU75N125 (multiplied by a factor 5), CusoZnso, and CU50AU50. 
Diamonds refer to the calculated values, the full line (back Fourier transform) serves 
as a guide to the eye 

exponentially damped. The different behavior of the amplitudes can be re- 
lated to differently large disorder in the neighborhood of relevant extremal 
points of the alloy Fermi surfaces. 

3.6 List of published applications 

We briefly review applications of the formalism developed in previous sec- 
tions to specific problems. Additional details concerning formalism and not 
discussed here in details, e.g., the expansion of the IEC expression in terms 
of the small parameter 1 — cos(£?) or the details concerning the numerical 
verification of the vertex-cancellation theorem, can be found in re- 
spectively. The influence of surface roughness (fluctuating spacer thickness 
and diffusion at the interface between spacer and magnetic layers) on the 
oscillation amplitudes was studied in Q . The effect of alloying in the spacer 
[ p"9| on the oscillation periods and their amplitudes, and in magnetic layers 
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[ pl| on the oscillation amplitudes and phases was also studied in detail for the 
trilayer system Co/Cu/Co(001). Ordering in disordered spacers ]22j and/or 
magnetic 1 20 layers lead to a formation of new periods not present in ideal 
spacers. Oscillations of the IEC can originate not only due to the spacer but 
also from adlayers or cap layers. We refer the reader interested in this prob- 
lem to a recent detailed study |23|,^4j . Finally, the study of the temperature 
dependence of the IEC and of the combined effect of the temperature and 
disorder is subject of very recent papers 0,|l7|], respectively. 



4 Conclusions 

We have derived closed expressions for the exchange coupling between two 
magnetic subsystems separated by a non-magnetic spacer with a relative an- 
gle 9 between the orientations of the magnetizations in the magnetic slabs. 
The derivation is based on a surface Green function formalism. The numerical 
effort scales linearly with the thickness of both the spacer and the magnetic 
slabs. The formulation allows also for an efficient evaluation of the tempera- 
ture dependence of the coupling amplitudes. Numerical examples were cho- 
sen to illustrate the theoretical aspects rather than to give a comprehensive 
overview of results obtained by the present formalism or by related methods. 

We wish now briefly to mention some unsolved problems. The following 
list is neither complete nor are the problems listed according to their im- 
portance: (i) The oscillatory dependence of the IEC on the thickness of the 
magnetic slabs was not yet systematically investigated on an ab initio level. 
Existing calculations |]l^ , ^8|p9tl were performed for too thin magnetic slabs to 
relate occurring oscillations to extremal points of spin-polarized Fermi sur- 
faces; (ii) The problem of biquadratic and higher order terms also did not 
receive a proper attention on an ab initio level. A relevant problem is a sys- 
tematic study of situations for which the non-collincar (biquadratic) coupling 
can dominate. Obviously, it can happen most probably for the spacer thick- 
nesses for which the IEC values are close to the transition between the F and 
AF couplings |fti[ . In addition, it remains to be seen whether a theoretical 
description of biquadratic coupling has to be based on a fully relativistic spin- 
polarized level; (iii) The stud y of superstructures in the spacer and/or in the 



magnetic slabs (see Sec. 2.10) offers a possibility of a deeper insight into the 
physical nature of the IEC because of new periods, which are connected with 
the extremal vectors of the spacer material in a more sophisticated manner 
than in the canonical cases of Cu or Cr spacers; (iv) The study of oscillatory 
behavior of exchange interaction across magnetic spacers is of great inter- 
est. One possibility here is to employ the method of infinitesimal rotations 
[||Q; (v) The study of exchange coupling through the semiconducting or, 
more generally, through an insulating spacer where one expects exponential 
rather than ./V~ 2 -decay has remained limited until now to model studies 
(vi) The study of alloying in the spacer, magnetic layers and at interfaces 
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has to be extended to new interesting systems. It offers a straightforward 
method to obtain valuable informations concerning alloy Fermi surfaces, in 
particular for the case of alloyed spacers; and, finally (vii) The study of the 
IEC through spacers with complex Fermi surfaces, in particular through the 
transition metal spacers. 
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A Vertex cancellation theorem 

We present here a general discussion of exchange interactions in the presence 
of substitutional disorder. The results given here are used in the present pa- 
per to study interlayer exchange interactions, but they are also applicable for 
studying exchange interactions within a ferromagnet, exchange stiffnesses, 
spin-wave energies, etc. The principal result is the "vertex cancellation the- 
orem" of Bruno et al. [|l5| . In here we give an alternative, more general, 
derivation of this result. 

Let u = {ur} be a particular configuration of the local moments, where 
ur is a unit vector pointing in the direction of the R-th local moment. We are 
interested in the variation of the thermodynamic grandcanonical potential 

1 f +0 ° 

J2 a = --Im/ f(E,T)Tr(ln gii (E + M + ))dE (60) 

^ J — oo 

with respect to u. The Green function ga(z) for a particular alloy configura- 
tion is defined from the potential function Pa(z) corresponding to u as 

9a(z) = (Pu(z) - Sy 1 . (61) 

An immediate consequence of (^lj) is a trivial commutator relation to be used 
below, namely 

[Pu(z);g a (z)}_ = [S;g a (z)]_ , (62) 

where [A; B]_ = AB — BA. The configuration averaged Green function 
(ga(z)} = 3u(^) is usually formulated in terms of the coherent potential 
function Vu(z) as 



Vu(z) = (Vu(z) ~ sy 1 



(63) 
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which leads to a relation analogous to (|6^) , 

[P a (z);g a (z)]_ = [£;&(*)]_ . (64) 

In general, the averaging in (^0|) cannot be reduced to In g^(z) and an evalu- 
ation of the so-called vertex corrections is necessary. We shall show, however, 
that the variation of ( |60| ) due to an infinitesimal change of u takes a simple 
form. 

Let us consider the variation of the potential functions Pa(z) in more 
detail. To each lattice site R we associate a non-random vector 0r = 6*r Ar, 
where Ar refers to the axis of rotation and #r to rotation angle by which the 
reference orientation u 0; r is transformed into Ur. The transformed potential 
functions are therefore given by the following similarity transformation 

where the rotation matrix /7© in (|6^) is defined as 
( U ®)nLs,WL<s> = 5 R,R' S l,l< x 

cos 1 ~ * sin " R (T ■ ( 66 ) 

The symbol er in (|66|) denotes the vector of the standard 2x2 Pauli matrices 
and 1 is the 2x2 unit matrix. The first-order change of Pa.n(z) caused by 
an additional infinitesimal rotation <5vr is then expressed as 

SP^z) = [(SA-;P a (z)]_ , (67) 

where the matrix elements of the operator SK = Us v — 1 are explicitly given 

by 

(^ri.,rw = d RR' 8l,l> ^ [a- ■ 5v R } SiS , . (68) 

The introduced infinitesimal rotation vectors (5vr satisfy Ug v U& = U&+s@ 
whereas, in general, Us& U@ ^ J7©+,5©. Let us note that 5K is a non-random 
site-diagonal operator. 

The first-order variation of Tr (In ga(z)) can be now formulated using ( p4| , 
|3) as 

STr (hi g a (z)) =-Tr(g a (z) [5K;P a (z)]_) , (69) 

which can be rewritten by applying the permutation invariance of the trace 
and (§|, H) as 

<5Tr(ln fffl (z)) = -Tr {SK ([P a (z)- g a (z)}_)} 

= -Tv{SK ([S;g a (z)]_)} 

= -Ti{5K [S;g^)L} 

= -Tr{SK[Vn(z);g^)}-} ■ (70) 
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By using the permutation invariance of the trace once again, (70) can be 
given the final form 



<5Tr(ln 3£l (z)) =-Tr{g a (z) {6K;V<x(z)]_} 



(71) 



Let us note that (71) was derived in a formally exact alloy theory, but is valid 
in the CPA as well. Within the CPA, the result ([n]) has an obvious interpre- 
tation: the r.h.s. describes the variation of Trlng a (z) induced by performing 
on the site-diagonal coherent potential functions Vu : r(z) the same rotations 
( |68| ) as applied to the potential functions -Pu,r(-z); note however, that this is 
not equal to the infinitesimal change of the true self- consistent CPA coherent 
potential function. 

Thus, the torque acting on the moment at site R due to the exchange 
interactions is given by 



u,R 



Sf2 a 
<5vr. 

Im Tri g Q (E + iO 



f(E,T) x 

. H) 



[n R <r; V a (E + iO+)] } dE , (72) 



where I1r is a projector on site R. This exact result constitutes the "vertex 
cancellation theorem" for the torque. Its usefulness arises from the fact that 
the "vertex corrections" have been eliminated. 

In order to compute the difference of thermodynamic grandcanonical po- 
tential between two local moment configurations Ui and U2 in the CPA, we 
use a theorem due to Ducastelle |31 , which states that the thermodynamic 



grandcanonical potential, considered as a functional fi\P,P] of the indepen- 
dent variables V and P, is stationary with respect to V when the latter 
satisfies the CPA self-consistency condition. This means that a first-order 
error in Va gives only a second-order error in i?u. Let us approximate Va(z) 

by 



(73) 



i.e., we assume that Va{z) is transformed like Pa(z) under a rotation of the 
local moment direction. This can be expected to be a good approximation, 
provided the condition 



d© 



R 



dR 



< k F q R 



(74) 



is satisfied, where q-R and raR are respectively the charge and spin moment 
at site R. We then get 



5 S W»5«W = (^)-S) • 
Replacing Va by V'^ and g Q by in (jn]), we obtain 

STt (]ngti(z)) » <5Tr \ng^(z) , 



(75) 
(76) 
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and integrating over the angles, we get 



+ OC 



f(E,T)x 



Im Tr [ln^ (E + i0+) - In (E + i0+)] dE , 



(77) 



which constitutes the "vertex cancellation theorem" for exchange energies. 
Note that we have derived here a form of the "vertex cancellation theorem" 
within the CPA since this is the scheme which is used in practical calculations; 
however, one can prove that the same result holds if one takes the exact 
solution to the configuration averaging problem. 

In the case of interlayer coupling, the condition (^) is satisfied even 
for large rotation angles, because d0R,/dR differs from zero only in a region 
where ttir is negligible. This was confirmed by explicit numerical calculations 
in H. 



B The interface-interface part of the grandcanonical 
potential 



In this Appendix we derive the basic relations for an evaluation of the IEC 
within the interface-interface interaction formulation. 

The subsystems L and 71 can be downfolded using the formula (p8[) 



Trln(P- S) =Tr £ ln P-S +Tr TC ln P-S 



+ Tr c In 



(P - S) CC ~(P- S)cc 



C 



P-S 



(P - S) CC 



(P-S) 



K 



CR 



P-S 



(P-S) 



RC 



(78) 



The first two terms are independent of the rotation angle 9 and, consequently, 
they do not contribute to the exchange energy £ x (9). We are thus left with 
a quantity which is limited to the subspace C only. It is now easy to identify 
the individual terms in (|7^) . The potential function blocks between different 
subspaces such as Pec or Pen are zero because the potential function P is 
site-diagonal. The blocks of 5* between neighboring subspaces do not vanish, 
but the non-zero subblocks connect only neighboring principal layers. The 
important part of the Tr In (P — S) is then reduced to 



Tr c In (P - S) CC + Tr c In 



I - 



P 



C C 
— ^SiqGcSqi — — — ^Soi^tcSio 



P-S 



(79) 



The first term is independent of 9 and thus does not contribute to the ex- 
change energy. The second term can be simplified using the two-potential for- 
mula @. We identify G (0) = C/(P-S), v x = SwGcSoi, and v 2 = S 01 g n S w - 
The t-matrices are then identical with the r-matrices, and the potentials v\ 
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and V2 are equal to the embedding potentials 7\ and ZV In this way we find 
the expression for the grandcanonical potential 



f(E,T,fi) Tri In 1 - gi N (z)T N (z)g N i(z)Ti(z) 



-— Im 

7T 



dz . 



(80) 



where f2o{T,iJ,) contains all the terms independent of 9 and the Tri applies 
only to the layer 1, i.e., the first spacer layer. If the system is invariant with 
respect to translations in the planes of atoms, or, if such a symmetry is 
restored by configuration averaging, (80) can be written as 



tt(0,T» = O)(T,/i)-ilm / f(E,T,n)x 

^2 trln 1 -5lA f ( k ||, 2 ) T A f ( k || , z )ffJVl(k||,2:)Ti(k||,z) 



dz , 



(81) 



where tr means the trace over angular momentum indices L = {£m) and the 
spin index a. 



C Useful mathematical tools 

Theoretical developments and many calculations are facilitated by the parti- 
tioning technique and the two-potential formula applied to the Green function 
and its logarithm. 

Let P and Q denote projection operators onto the complementary sub- 
spaces (i.e. P + Q = 1). We denote the projections of matrices as PAP = 
App, PAQ = Apq, etc., and P/A means the inversion of App in the subspace 
referring to projector P. In most applications, A — z — H ot A — P{z) — S 
and G{z)=A-\ 

The projections of the inverse A -1 to the matrix A are given by |Q 
{A " )PP = A A *A ' (82) 

App - Apq-^Aqp 

(A-') QQ = - Q P , (83) 

A QQ - Aqp^Apq 

(A- v )p Q = -jA PQ {A- 1 ) QQ = -{A-^ppApqQ , (84) 
(A-'hp = -^A qp {A- 1 )pp = -{A-^qqAqp^ . (85) 

It is sometimes easier to invert the full matrix A than its blocks. In such 
a case the inverse partitioning is useful 

r (A-I)pp - - ^ = 5 A (A-I)pp . (86) 



Ipp ' P - A PQ {A-i) QP P-{A-i)p Q A QP 
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This can be used to calculate the surface Green function of a semi-infinite 
system from the Green function of the infinite system. 

Partitioning technique also allows to simplify calculations involving Tr In 
of a matrix. The basic relation is 

Tr In A = In det A. (87) 

It then follows Tr In AB = TrlnA + TrlnS, Trlnl = 0, Trln (A -1 ) = 
-Trln A, and Trln [(A - B)^ 1 ] = -TrlnA - Trln[l - A~ 1 B}. The Tr In 
A can then be partitioned as 

P 

Tr In A = Tr P In [PAP] + Ttq In [QAQ - QA—AQ] . (88) 

J\ 

To prove (|88|), let us multiply the matrix A by L = 1 — Aqp(P/A) from left 
and by R = 1 - (P/A)A PQ from right. The result is LAR = A PP + A QQ - 
Aq P (P/A)A P q. Now using @, and the fact that det [L] = det [R] = lwe 
find (jS8|). In a special, but important case, when A PP = P and Aqq — Q it 
holds 

Tr In A = Tr P+Q In [P + Q + A PQ + A QP ] 

= Tr P In [P - A PQ A QP ] = Tr Q In [Q - A QP A PQ ] . (89) 

The Green function of a system described by the Hamiltonian H = H + 
vi + v 2 , where H is the unperturbed part, and Vi(i = 1,2) are perturbing 
potentials, is given by G = G<® + G^TG (0 \ where G = (z — H)' 1 , G<® = 
(z - Hoy 1 , and T = V(l — G^V)~ X , where V = v t + v 2 . The full T-matrix 
T can be expressed in terms of the t-matrices, ti — Vi(l — G' '^) -1 , (i = 1, 2) 
and of the unperturbed resolvent G^ by the the two-potential formula 

T = ti[l- G^hG^h]- 1 (1 + G^t 2 ) +t 2 [1 - G^hG^h]- 1 x 

(1 + G (0) ti). (90) 

It is derived in the following way. Because 

(1 - A)[l - (1 - A)" 1 AB(1 - B) -1 ](l - B) = 1 - A - B , (91) 

it holds 

Tr In [1 — A — B] = Tr In [1 - A] + Tr In [1 - B] 

+Trln[l - (1 - A)- 1 AB (1 - B)' 1 } . (92) 

By inserting A = G^v% and B = G^v 2 into @ one obtains @. The 
two-potential formula for the Tr In of the full Green function 

TrlnG = TrlnG (0) [l - VG^]' 1 

= TrlnG (0) - Trln [1 - G ( %i - G (0) v 2 ] 

= TrlnG (0) - Trln [1 - G ( %i] - Trln [1 - G {a) v 2 ] 

-Trln[l-G (0) tiG (0) t 2 ] (93) 
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follows directly from (|92[). 

If the matrix A is a function of a variable z (complex in the general case), 
the derivative with respect to z is given by 



provided that the matrix A(z) is nonsingular. This identity is used to derive 
the expression of the grandcanonical potential fi in terms of the auxiliary 
Green function @ within the TB-LMTO. 

The identity in (|87]) is valid up to an integer multiple of 2m. Neglect- 
ing this fact can lead to serious errors. There is no panacea for this kind 
of difficulties, but in some situations they can be avoided, for example by 
choosing the integration contour parallel to the imaginary axis, but this is 
not always possible. In some cases the incremental procedure for calculating 
the In det, \a f{zk+i) = ln/(zfe)+ln [f (zk+i) / f (zk)} in the spirit of an analyt- 
ical continuation can be helpful, provided that the change of phase between 
two consecutive points Zk is less than 2-7T. To insure this, one has to choose a 
sufficiently small grid in z. 

D Inversion of block-tridiagonal matrices 

We wish to compute g = A^ 1 for a block-tridiagonal A. The matrix A 
is divided into N x N square subblocks of the same dimension m, from 
which non-zero are only Ak,k, ^4fe-i,fe, and Ak^k-i- The diagonal blocks are 
a sum of two terms: hermitean matrix and a symmetric complex matrix. 
They are always non-singular. The off-diagonal blocks under the diagonal 
are equal to hermitean conjugate of the corresponding blocks above the diag- 
onal (Ak y k-i = -^t-i &)• The methods based on repeated use of partitioning 
are particularly efficient if only diagonal blocks, or four so-called 'corner' 
blocks (51,1, gN.N, Si, ./v, 9n,i) are needed like in the interlayer exchange cou- 
pling calculations. 

First, four sequences of auxiliary matrices are calculated 

Xn-U — AN-k,N-k+l(AN-k+l,N-k+l — Xjy^k+l)~ 1 AN-k+l,N-k , 



-^-TrlnL4(z)] = Tr \-^A(z) A _1 (z) 



(94) 



Y k+ i 



X N = 0, (k = l,...N-l), 

Ak+i,k(Ak,k — y'k)~ 1 AkM+i , Yi 
-(A^k - Xk)- 1 A^k-i, (fc = 2, 
-{Ak.k-Yk^A^k+u (k = l,. 



0, (k = 2, 
■AO 

N-l), 



N) 



(95) 



that are used to compute the diagonal and off-diagonal blocks of g 

gk.k — {Ak,k — Xk — Yk) 1 , 
9i,j = g, : . , for i > j , 
9i,i = Wi gi+i,j for i < j . 



(96) 
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It can be proved that the numerical effort to evaluate the corner blocks scales 
as 0(Nm 3 ). The details, particularly the tests of efficiency can be found in 
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